Moment distributions of clusters and molecules in the adiabatic rotor model 
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We present a Fortran program to compute the distribution of dipole moments of free particles for 
use in analyzing molecular beams experiments that measure moments by deflection in an inhomo- 
geneous field. The theory is the same for magnetic and electric dipole moments, and is based on a 
thermal ensemble of classical particles that are free to rotate and that have moment vectors aligned 
along a principal axis of rotation. The theory has two parameters, the ratio of the magnetic (or 
electric) dipole energy to the thermal energy, and the ratio of moments of inertia of the rotor. 
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I. INTRODUCTION 

It is common to measure the magnetic moment or 
the electric dipole moment of clusters or small molecules 
by the deflection of a molecular beam by an inhomoge- 
neous ficld[l-5]. These experiments take place in the gas 
phase using Stern-Gerlach magnets to deflect the beam 
in the magnetic case[4-5] and using an electric field gra- 
dient in the electric case [1-3]. One needs a theory of 
the moment distribution to relate the observed deflec- 
tions to the intrinsic moment of the particles. There 
are two well-known limits for the distribution, the quan- 
tum limit of a spin with a fixed total angular momen- 
tum, and superparamagnetic limit, where the moments 
are thermally distributed. Neither of these limits is valid 
for the typical situation of a nanoparticle, which may 
have a moment with fixed orientation in a body-centered 
frame, but changing orientation in the laboratory frame. 
If we assume that an external field is introduced adia- 
batically, the distribution of moments can be computed 
using the adiabatic invariants of the rigid rotor. This 
classical theory and a method of solution was given in 
ref. @. While conceptually the theory is quite straight- 
forward, the computation is not completely trivial. Wc 
first present the equations that govern the deflections, 
and then the computational aspects. 



II. THEORY 

We consider a ferromagnetic particle having a mag- 
netic moment aligned along the 3-axis and equal mo- 
ments of inertia around the 1 and 2 axes in the body-fixed 
frame. Its Lagrangian is given by 



L= — 



i> 2 sin 2 I 



+ /jqB cos 9 
(1) 

where J\, J 2 = Ji and J3 are the principal moments of 
inertia and 0, <f>, ip are the Eulerian angles of the 3-axis 
with respect to the magnetic field B. The theory would 
be the same for a particle with an intrinsic electric dipole 
moment po in an electric field E, simply replacing [i^B 



by poE in all equations. There are three constants of 
motion for the Lagrangian eq. (1). They are the energy 
E, 



b 2 sin 2 I 



2 



(ip + cf> cos 9^j —(IqBcosO 
the angular momentum about the field direction m z , 



dL 

— - = <A< 

d(j) 



+ J3 ( ip + i> cos & ) cos 



and the angular momentum about the 3-axis 7713, 



dL r - 

777.3 = r = J3 ( 1p + (p COS ( 

dip 



(3) 



(4) 



The last quantity, 7713, is only conserved because of the 
condition we imposed that J2 = J\. Under the equa- 
tion of motion, the variable 9 has a periodic dependence 
on time, oscillating between two limits 9\ and #2- For 
convenience below, we replace the variable 9 by its co- 
sine, u = cos#. The quantity of interest for the deflec- 
tion measurement is the average moment of the particle 
/2 = noii where u is the average of u over a cycle. There 
is an analytic formula for this quantity in terms of the 
elliptic integrals K{y) and E(v) which can be compactly 
expressed in term of the three zeros Uo, U\, U2 of the cubic 
polynomial 

f(u) = (2Ji-E-Jim^/J 3 +2Ji/i 5u)(l — Ti 2 ) — (m z — 777,377) 2 

_ (5) 
The formula for u is 1111] 



u K(v) + (u 2 - u )E(v) 
K(v) 



(6) 



where v = (u 2 — ui)/(u 2 — Uq). 

We use eq. ((HJ) to compute u as a function of 7773, m z 
and E. However, E changes as the particle enters the 
field. Assuming the field change is adiabatic, the action 
Jg associated with the variable 9 remains constant and 
thus can be used to determine the new value of E. There 
is no analytic expression for E{Jq) or even for the inverse 



2 



function Jg(E). In the program we compute the latter 
from its definition 



1 - u 2 



(7) 



In zero external field, the action is simply related to the 
total angular momentum J, 



I = max(|m 3 |, \m z \) + 

Z7T 



(8) 



This relation is useful to make a connection to the quan- 
tum mechanical formulation of the problem as well as to 
make tests of the program. 

The probability distribution P(u) that we seek to com- 
pute can now be expressed as the three-dimensional in- 
tegral, 

P(u) = 



Z{T) 



dl 



-i 



dm z / dm,3 5(u—u(I,m Zl m3))e E °/ kT 



-/ 



(9) 

where the partition function Z(T) is the corresponding 
integral without the delta function and Eq = (I 2 — 
m 2 )/2Ji + 777.3/2 J3. There are two symmetries that 
can be used to reduce the size of the integration re- 
gion. Namely, 77(7,7772,7773) remains the same under the 
interchange of m z and 7773 and under the replacement 
ra3,m 2 — > — 7773, — m z . While eq. ([9]) is expressed in 
terms of dimensioned physical parameters, in fact the re- 
sults only depend on two dimcnsionlcss combinations of 
those parameters, namely 



kT 



(10) 



and J\l J3 Note that the distribution function is indepen- 
dent of the overall magnitudes of the moments of inertia. 



III. NUMERICAL 

We evaluate the integral ^ using uniform meshes in 
the three integration variables, binning values of 77 on the 
mesh points to construct the probability density. This 
requires a fine integration mesh due to the singularities 
and discontinuities in the integrand. We use a mesh size 
of Am/ 1 « 0.005 to achieve an accuracy suitable for 
graphing the distribution P{u). It also helps to have 
incommensurate mesh spacings for two m integrations. 

Another numerical problem is connected with deter- 
mining u as a function of Jg. Both quantities are com- 
puted directly in terms of the energy variable E, but to 
find 77 as a function of Jg requires solving an implicit 
equation. In the program this is carried out by Newton's 
method; a warning is given if the convergence is poor. 




FIG. 1: Zero field distribution compared with the numerical 
results for x = 0.01 and J1/J3 = 1. 



IV. TESTS AND PROGRAM USE 

There are two analytic tests that can be made of the 
program. The first is the probability distribution at zero 
field, which is given [8[ by 



p(«) = iiog(i/M). 



(ii) 



Unfortunately, eq. (6) can not be used at B = because 
uq goes to infinity at that point. However, the numerical 
parameters in the program have been set so that the dis- 
tributions are accurate to within a few percent for values 
of x greate than 0.01. Fig. 1 show the comparison of 
eq. (|11|) with the computed distribution at x — 0.1 with 
the mesh as given above. The small irregularities are the 
binning effects associated with the finite mesh size. 

The second analytic test is the ensemble-average mo- 
ment (77) at small fields. It is given by 



(u) 



(12) 



The computed ensemble average for x = 0.01 is (u) — 
0.000222, in excellent agreement with eq. (|12jl . 

The program runs without any input file, as all of the 
parameters have been set in the Fortran coding. The 
important physical parameters x and J\ / J3 are specified 
on lines 22 and 25 of the code, respective. Running the 
code with the values given, 

betamuOB=l . OdO 
JlJ3=l.dO 
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gives as direct output the values of x, J3/ J\ and (it), 



1.000 1.000 0.19220 



The program also writes a data file 'udist.dat' that has 
a table of values of u and P(u). Fig. 2 shows a plot of 
that data. 



FIG. 2: Distribution P(u) computed for x = 1 and J1/J3 = 1. 
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